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Abstract. To ensure homogeneous conditions within the complete area of perfused 
microfluidic bio- reactors, we develop a general design of a continuously feed bio-reactor 
with uniform perfusion flow. This is achieved by introducing a specific type of perfusion 
inlet to the reaction area. The geometry of these inlets are found using the methods 
of topology optimization and shape optimization. The results are compared with two 
different analytic models, from which a general parametric description of the design 
is obtained and tested numerically. Such a parametric description will generally be 
beneficial for the design of a broad range of microfluidic bioreactors used for e.g. cell 
culturing and analysis, and in feeding bio-arrays. 
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1. Introduction 

The development of microfluidics, to handle minute amounts of fluids, is currently 
revolutionizing fluid transport in the field of analytic cell-biology: Traditionally, cells 
are cultured in so-called batch cultures in a flask and an experiment is typically initiated 
by adding an agent. After a certain time, such as a day or two, the response of the agent 
is studied using typically only one reporter such as fluorescence. In order to increase 
throughput, cells can, at present, be cultured and assayed in robotically controlled 96 
or 384 well plates. By contrast, culturing of cells on a microfluidics device gives a range 
of new possibilities [1] e.g. studying cell mobility in real time when exposed to stable 
continuous gradients [2j. Furthermore, combinatory experiments can be performed on 
chip that are based on arrays of interconnecting chambers [3l H] . 

The inlet design presented in this paper introduce a number of improvements to 
current perfused bio-reactors: The creation of uniform flow conditions all over the bio- 
reactor ensures homogeneous cell conditions both with regards to concentrations of 
externally supplied growth-factors and to the shear induced on the cells by the perfusion 
flow. Too small a height of cell culture chips is inhibiting cell growth [5], El [7] , and in 
Refs. [HI [9] it has been shown that the chamber height must exceed 1.5 mm in order to 
provide identical culturing conditions as in traditional cell culture flask. On the other 
hand, to ensure laminar flow conditions, a small height is preferred. Therefore in the 
case of cell-culturing chips, a chamber height of 1.5 mm is optimal. In other cases, 
such as for micro-array hybridization chambers, the functionality indeed benefit from 
far smaller reactor-channel heights, where volume needs to be minimized in combination 
with a maximization of reaction area. 

In recent years different bio-reactors have been constructed, where the uniformity 
of the perfusion flow along the reaction area has been achieved at the expense of a large 
hydraulic resistance across the whole bio-reactor [3]. One example is the Micro cell- 
culture chamber by M. Stangegaard et al. [8], where the fluid is directed from a wide 
reservoir through a large number of small parallel channels. This barrier creates a large 
pressure drop which give rise to the uniform flow. From the inlet structure described 
from this work, the same uniform flow-field can be achieved with a significantly lower 
pressure- drop. This opens up the possibilities of driving the perfusion flow by low- 
power methods such as e.g. buoyancy force, which recently has been used to drive other 
microfluidic devices [18j. 

Our novel design also reduces the fluid volume used in creating the uniform flow, 
which is crucial both when dealing with expensive biochemical samples or to avoid 
dilution of small samples. Additionally it will enable a better analysis of fast cell- 
reaction kinetics with high time resolution. 

The paper is organized as follows: In section [2] the general bio-reactor layout is 
outlined together with an introduction of the related characteristic parameters. In 
section [3] the optimal structure of the perfusion inlet is found, first by the general 
method of topology optimization, which imposes no constraints on the topology of the 
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Figure 1. A sketch of the bio-reactor with its three sections: (I) the perfusion inlet 
of constant height hi, (E) the expansion chamber of spatially varying height he(x,y), 
and (II) the main reactor area of constant height /i2- The perfusion flow, driven by 
the pressure drop Ap, is indicated by arrows. 



structure. The resulting structure is further refined by the method of shape optimization. 
The optimized geometry is in section H] compared to two simple expansion design, while 
in section [5] the results are summarized in a design guide. The analysis and the design 
guide is further verified by full 3D simulations in section [61 A conclusion is given in 
sec. 

2. Geometry and basic flow equations of the microfiuidic bio-reactor 

The generic microfiuidic bio- reactor layout used in this work is illustrated in figure [H It 
consists of a single microchannel perfusion inlet (I) of constant height hi, which broaden 
out in an expansion chamber (E) of varying height he{x,y) to distribute the fiuid over 
the much wider and more shallow main reactor (II) of constant height /i2, where the 
cells are immobilized. All vertical channel and chamberheights in the z direction are 
much smaller than any lateral length scale in the xy plane; the bio-reactor is thus fiat. 

The main objective is to obtain a uniform fiow in the main reactor with minimal 
pressure drop Ap and with a minimal volume of the expansion chamber. This is achieved 
by carefully designing variations in chamber height he{x,y) of the expansion chamber. 
As the constant inlet channel height hi is assumed larger than the constant height of 
the main reactor area /12, the height variation in the expansion chamber (E) will be 
bounded by these two heights: 



The whole bio-reactor is assumed symmetric both through a central vertical and a 
central horizontal axis, and as a consequence only the upper left part will be dealt with 
here. 

As we consider only low concentrations of the solutes and a constant temperature, 
the density p and viscosity t] of the buffer liquid are constant in space, and the fiow is 
determined by the geometry of the reactor and the applied pressure drop Ap driving the 
fiow. As a consequence of the assumed fiatness of the bio-reactor, the pressure p does 
not vary in the vertical z direction, i.e. dp/dz = 0. Moreover, due to the small heights, 
viscous damping from the top and bottom plates of the bio-reactor dominates the fiuid 




(1) 
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flow and makes the flow laminar. This is evident from the value of the Reynolds number 
Re given the low flow velocities, u ^ 1 mm/s, and small length-scales, h ^ 1 mm of the 
system: Re ~ 1. 

In this flow regime it is useful to work with the 2;-averaged 2D velocity field u{x, y) = 
{l/h{x,y)) w{x,y, z) dz of the full 3D field v(a;, y,^;). To a good approximation u 
fulfills the 2D Brinkman-Darcy equation [TO] . 



12n 

V'u - j^^^^ u - Vp{x, y) = 0. (2) 
Here the prefactor 12'r]/h'^{x,y), also denoted the damping coefficient a, 

is reminiscent of the z-part of the Laplace operator in the full 3D-description, and it 
represents the dominant part of the viscous damping of the liquid in the system. 

For possible continues changes in the height h{x,y) of the expansion region (E), 
the z-averaged 2D velocity field is not divergence-free due to mass-conservation, but 
an additional tern arises: V ■ u{x,y) = —\'Vh{x,y)\/h{x,y). In the case of possible 
discontinues jumps in height along an interface, this correction becomes the following 
new boundary condition on the interface: 

hiUi ■ n = h2U2 • n, Ui • t = U2 ■ t, pi = p2, (4) 

where the height-subscript is extended to the corresponding velocities and pressures. 

Working with the 2D-restricted description, the detailed geometry of the expansion 
chamber is illustrated in figure [2l The important in-plane length scales are the length L 
of the expansion chamber, the width W of the main reactor, and the width ^ of the inlet. 
To achieve a uniform flow in the main chamber, the pressure along the line A dividing 
the expansion chamber and the main chamber must by constant. Consequently, the 
spatial variation of the expansion chamber height he{x,y) must be optimized in order 
to get as homogeneous pressure along line A as possible. 



3. Optimization 

To enable optimization of the system two additional parts are introduced. First, a 
set of design variables 7, which uniquely characterizes all available configurations in 
the optimization problem, and for which a unique solution U(7) to the system exists. 
Second, an objective function $ which quantifies how well a given configuration of the 
system performs. By convention $ has to be minimized in order to achieve the optimal 
solution, and generally the objective function can depend on the design variables and 
the related solution of the system $ = $(7, U(7)^. 

As alluded to in the previous section, we base our objective function on the 
homogeneity of the pressure along the cross-section A, since a uniform pressure there 
will lead to the required uniform flow field in the main reactor. In the following, this 
objective will be expressed in two different ways, depending on the given optimization 
methods. 
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Figure 2. (a) The geometry of the expansion chamber E. The parameters related 
to geometry and flow are defined together with specific segments and points used for 
later reference. The optimal transition between the inlet height hi (white area) and 
the main reactor channel height (gray) occurs inside the expansion region E (light 
gray), (b) Topology optimized height variation in the expansion region. The height 
changes almost step-like from the white region F (hi) to the gray region D (/i2)- Note 
the small region of intermediate height marked G. 



3.1. Topology optimization of the spatial height variation 

To search for the globally optimal solution, and not a priori exclude any non-intuitive 
solutions, we will not rely on any pre-described variation of the height. Therefore, 
we begin by applying the method topology optimization [12], which by definition is 
independent of the topology and therefore unlimited in its search for the optimal bio- 
reactor design. The method of topology optimization was first applied to the field of 
structural mechanics [13], and have been recently implemented to the field of microfiuidic 
systems [TT| [H] and chemical microreactors [15] . 

Arbitrary height variations of he{x, y) can be realized by representing the height as 
a variation of the design variable field '~f{x,y), where < 7 < 1. To cover the range 
of heights defined in equation ([1]), the design variable is assigned the value 7 = to 
describe chamber heights equal to the inlet height hi and the value 7 = 1 for heights 
equal to the main reactor height h2. In the expansion chamber, now denoted the design 
region Q, the design field can take any value < 'y{x,y) < 1 to describe all possible 
height variations he{x,y) . 

The actual implementation, method and procedure of topology optimization will 
not be touched upon here, as it is fully described in the work of Olesen, Okkels and 
Bruus [Hj. Still what is essential for this work is the objective function $, which has 
to be chosen with care. To obtain a numerically stable search we define $ as the square 
deviation of the pressure around a reference pressure Pref along A: 

^ = ^ Ijp - p,,,)Ms. (5) 

We choose pref as the pressure at the far corner C in the case where the expansion 
chamber has the same height hi everywhere as the inlet. 

Figure [2]^b) shows the resulting optimal height distribution hf>{x, y) for the following 
set of parameters: W = lO'^ m, L = 0.95 W,i = 0.1 W, hi = 0.04 W, and /12 = 0.02 W, 
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where the gray-scale color-coding spans from height hi in white down to /i2 in gray. We 
define the ratio A between the two damping coefficients 12ri/h\ and 12ri/h\ as 



and get the value of ^ = 4 for figure [2](b). 

As mentioned earlier, any change in height produces a correction-term to the 
continuity equation, when using a 2D-restricted description. It turns out that both 
ways of implementing this correction in the topology optimization problem fails, due 
to the very nature of the method. First, the free variations of the design field in 
topology optimization prohibits any interface to be defined a priori, and therefore the 
boundary conditions of equation (jlj) cannot be applied in this step of the optimization 
procedure. Second, it turns out that the solutions of the topology optimization problem 
involves sharp transitions in the height, limited by the grid-meshing length-scale of the 
finite element method. Therefore when including the correction-tern to the continuity 
equation, a fiuid source is added to single mesh-elements, and this destabilizes the 
convergence of the method. The way to work about this limitation, is to add the 
boundary conditions of equation (jl]) to the shape-optimization method, applied later in 
the optimization process, after the shape-optimization has been preliminarily compared 
to the topology-optimization. 

From the topology optimized solution in figure EJ^b) we see that among all possible 
height variations, the optimal design consists of a single sharp transition between a 
region of inlet height /ii, and a region of main reactor height hi. Only very close 
to the inlet channel is seen an ambiguity which indicate the possible existence of a 
region of intermediate height. From topology optimizations for other parameter- values, 
similar solutions arise with a sharp transition between regions of height hi and /i2; and 
consequently we conclude that such a single-connected transition is indeed the optimal 
solutions of the problem. 

To evaluate the quality of the topology optimized solution, we plot in figure [3] 
the pressure contour lines of the solution in figure [2](b), including a contour line 
corresponding to the value p = pref which goes through the corner C of the expansion 
chamber. Except close to the upper side wall, the pressure is seen to be uniform at the 
entrance of the main reactor, and it decreases uniformly throughout the whole extension 
of the main reactor. 

The chosen parameters used for the solution in figures [2](b) and [3] represent a rather 
extreme case, i.e., a combination of a small height-difference A = 4 and a wide expansion 
L/W = 0.95. As a result the transition extends nearly through the whole expansion 
region, but for all common parameters, the type of solution remains optimal. 

From the results of the topology optimization it is therefore natural to proceed with 
the shape optimization method, which compared to topology optimization involves fewer 
design parameters, is faster, and is numerically more stable. 




(6) 
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Figure 3. The contour lines of the pressure p{x, y) in the topology optimized design. 
The contour p — Prof , going through the corner C is marked. 



3.2. Shape- optimization 

In shape optimization the interface hne between the heights hi and /i2 in the transition 
chamber is given by a cubic interpolation hne through a number of control points 
{xi,yi),i = l,2,...,Nf3 as shown in figure a) with iV^ = 6. It is convenient to 
parameterize the points by the expression 

i 

{xi,yi) = {s,L, Si(3iW), — < Si < 1, < A < 1, (7) 

which ensures that all points lie within the expansion chamber. By fixing the factors jSi 
by Pi = [1 — 1/(3 iV/j)] X [{i — l)/[Np — 1)], a relatively even distribution of the control 
points is also ensured. During the optimization process the position of the interface line 
is changed by adjusting the control points Sj. 

The optimization is carried out by a simplex-method relying only on values of 
the objective function $(7), and not its partial derivatives d^/d^. To ensure efficient 
convergence of the given simplex method, it is beneficial to assign initial values around 
unity for the design variables 7. Furthermore, the method is unbound i.e. the design 
variables must give rise to a well-defined geometry regardless of their value. All this is 
accomplished by using the arcus tangent function: 

f 2 

s{(3i) = 1 - 7o (1 - A) 1 1 + - arctan 

with 

7 = {7o,7i,---,7iv^}- (9) 

We use Nj3 + 1 design variables to determine Np shape parameters because a faster 
convergence is achieved by adjusting the extend of the whole interface by a single 
parameter 70. Furthermore, equation ([8]) let the initial configuration of 

7i,it = [7o,l,...,l] (10) 

give rise to a well-defined, straight interface line reaching from the position (x, y) = 
(^{1 — 7o)-^v, 0^ to the upper corner. 



vr 



(7. - 1) 



z = l...Np, (8) 
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a 




s = ■ S = s(/33) S = 1 \ / 7v I 

Figure 4. (a) Setup for the parametrization of the shape-optimization problem. 
Variables s and /3 parameterize the normalized x and y-axis respectively, and the 
control-points s{Pi) are shown by circles with the interpolated interface curve in solid, 
(b) The pressure contours (thin lines) of the shape optimized positioning of the interface 
line (s, s/3) (thick line). Similarly to figure[3]the pressure contour are originating from 
the upper corner C 



Now that the interface by definition extends to the upper corner, this constraint 
does not need to be included in the objective function $ of the shape-optimization, and 
$ can therefore be defined with the sole purpose of achieving a uniform pressure along 
segment A: 



The actual optimization of the design-variables follows two steps: First a rough 
initial interface is found by using the initial setup of equation (ITO!) . and only adjusting 
the single variable 70, using a simple Matlab implementation of a bounded golden 
section search with combined parabolic interpolation [19]. Once a suitable straight 
initial interface is found, the actual shape is obtained using a direct unbounded simplex 
search method, also implemented in Matlab |20j . 

3.3. Results 

First, the shape-optimization method has to be validated with respect to the topology 
optimized solution, shown in figures [2](b) and[3l The same parameter- values were used, 
and the resulting shape-optimization shown in figure IH^b), is indeed similar to figure [31 
When comparing the results of the two optimization methods, it is observed that both 
the shape of the interface and the corner-pressure contour matches very well. Thereby 
we conclude that the shape-optimization is appropriate for the further analysis of the 
optimized interface. 

Again it should be noted that the set of parameters used in figures |2]^b), [3l and|4](b) 
is an extreme case, and therefore the shape-optimization method has been further tested 
to ensure the validity of this simple type of solutions. 




(11) 
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Figure 5. (a) A study of alternative expansion geometries. The expansion region in 
the Empty design is completely filled with inlet height hi, while the Box design has a 
simple box of height hi as expansion region, (b) Plot of the standard deviation Sp{s) 
of the pressure vertically across the expansion region as a function of the normalized 
horizontal position s into the region, as seen in (a). The different curves are topology 
optimized (solid line), shape optimized (dashed line), empty design (dotted line), and 
the four types of box design (dot-dashed lines), where the arrow mark the order of the 
values L = 0MW,0.55W,0.75W,0.95W. 



4. Comparison with alternative expansion geometries 

Now that the design of the expansion region has been optimized, it is natural to compare 
its efficiency to other alternative expansion designs. The ffist obvious candidate is to 
uniformly fill the existing expansion region with height hi i.e. to remove the topology 
optimization distribution of height h2, and this we will call the empty design. The next 
design comes as we replace the expansion region with a simple box of width W and 
height hi, and this will be called the box design. Both alternative designs are shown in 
figure [5](a). To compare these new candidate designs to the optimized designs, we will 
measure the homogeneity of the pressure around the end of the expansion region. To 
get an quantitative measure of the homogeneity of the pressure in the first part of the 
reactor, we measure the standard deviation 6p{x) of pressure across the width along the 
y-axis of the reactor part for a fixed x-coordinate: 

Sp{x) = ^{[p{x,y) - {p{x,y))y]^)y (12) 

where ( ■ )y is the mean along the y direction. Generally 6p{x) will decrease exponentially 
with the distance into the uniform reactor-part, and therefore 6p{x) will appear as an 
approximately straight line when shown shown in a log-linear plot as a function of x, 
see figure Mjo)- 

From the measurements presented in figure [5][b) it is first noted that both optimized 
designs produce a more homogeneous pressure-field than the alternative designs. While 
the empty design give the poorest results, the box design comes closer to the shape 
optimized design, and this tendency strengthen when moving the interface closer to 
the inlet e.g. for L = 0.35W. Since the box design evens out the pressure due to the 
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Table 1. The volume of the different designs of the expansion chamber, using the 
abbreviations: TO = Topology Optimized, SO = Shape Optimized, ED = Empty 
Design, and BD-X = Box Design, with the corresponding length-to-width fraction 
X — L/W. Second row shows the volumes in relation to the Shape optimized design. 



Type 


TO 


SO 


ED 


BD-0.95 


BD-0.75 


BD-0.55 


BD-0.35 


Vol (/iL) 


27.7 


27.6 


37.6 


68.4 


62.0 


55.6 


49.2 


Vol/Vol(SO) 


1.01 


1 


1.36 


2.48 


2.25 


2.01 


1.78 



translation invariant properties in the reactor part, there is a limit in how fast this can 
happen, as reflected in the slope of the dash-dotted lines in figure [5](b). On the contrary, 
the optimized designs aims at homogenizing the pressure by designing the expansion- 
parts, and therefore their corresponding slopes are steeper than the box design. As 
a result, the optimized designs are most efficient in quickly producing homogeneous 
pressure- fields. 

The hydraulic resistance i?hyd = ^p/Q of the expansion-regions for all the presented 
designs are in the range i?hyd = (1-7 — 3.7) x 10^ kgm~*^ s~^, which is five orders of 
magnitude smaller than the numerically estimated -Rhyd ~ 2.2 x 10^°kgm~^s~^ for the 
micro cell culture [8]. 

The optimized designs possesses another advantage, since the fluid-volume of 
the corresponding expansions regions are significantly smaller than any of the other 
designs mentioned. The fluid-volume of the different expansion-regions has been 
calculated/measured and is presented in Table [H Also presented in the table is the 
volumes relative to the Shape optimized design, and this clearly shows that extra fluid- 
volume is significantly higher especially for the box designs. 

From the above results we conclude that the optimized designs are generally better 
than the alternative designs, and we will therefore in the following present a general 
description based on a vast range of different shape optimized designs. 

Knowing now the basic shape of the height interface in the expansion region, we 
can now apply the right mass-conserving boundary conditions of equation (jlj) to the 
interface, and thereby improve the model upon which the following design guide is 
based. 

5. Design guide 

It is possible to match the numerically optimized geometry by simple theoretical models, 
which only depend significantly on the parameters W, L, A, as the remaining parameters 
i and hi only introduce minor corrections. By fitting the resulting interface obtained in 
these models for given parameters to the corresponding shaped optimized interface, we 
obtain an approximate parametrization which can serve as an easily applicable guide 
for practical design purposes. 

The basic idea behind the simple models is sketched in figure [61 Given the laminar 
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nature of the flow, we consider an idealized narrow flow stream stretching from the inlet, 
across the expansion chamber, to the entrance of the main reactor. The first part of the 
stream, which is the lower part of figure Et, goes along the horizontal symmetry- axis and 
starts at the inlet where the hydraulic damping factor is ai given by the height hi, see 
equation ([3]). Then, at the point {xo,y{xo)) the stream hits the interface, and continues 
horizontally to the point (L, y{xQ) with hydraulic damping factor 02- Along all streams, 
the hydraulic resistance is proportional to the effective length Leg = aiLi + 02/^2; where 



Li = yxQ + y{xQY and L2 = L — xq is the length of the first and second part of the 
stream, respectively. Since we seek the shape (xq, y(xo)) of the interface giving rise to the 
same pressure drop along the streamlines, all streamlines must have the same effective 
length. The specific form of the effective length, with its squares of Xq and y{xo), then 
leads us to expect an expression for y{xo), or in dimensionless form, an expression for 
s{(3) of the form 



The proposed expression for s{P) is of course not exact, but by calculating the shape 
optimized interface for a large number of parameter values, we can fit equation (fT3ll and 
make an statistical analysis of the obtained fitting parameters Sq^p and Cp. The resulting 
explicit parametrization becomes 



This parametrization is deduced for a plug flow model, which is shown in figure [6]d, 
while a more refined radial model, seen in figure Eb, can only be solved numerically. A 
comparison between the two models, is seen in figure [3 Here, the position s(0) of the 
interface at the center axis in the simple models for a large number of parameter values 
is compared to that of the shape optimized model. These results show no improvement 





(13) 




(14) 



(a) 




(b) 



(c) 



Li/W s{0)L L 



M/W s{0)L L 



Figure 6. (a) The concept behind the models where the flow along two distinct flow 
stream in the expansion chamber are compared, (b) The plug flow model, and (c) the 
radial flow model. Also shown are the involved parameters and constraints. The gray 
area corresponds to regions of height /12 • 
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0.8 0.85 0.9 0.95 

Theoretical sp(0), sr{0) 



Figure 7. A comparison of the calculated beginning s(0) of the transition region at the 
center axis between the simple model predictions and the shape optimization model. 
The comparison includes the simple model (crosses) and the radial model (circles) , and 
a perfect match would lie at the diagonal (solid line). 

by the radial model, and therefore it is adequate to base the design guide of equation 
[H] on the plug model. 

The design guide does generally a very good job, but it maybe worth adjusting 
the parameters used above (5o,p = 2.63, and Cp = 0.27) if the flow-homogeneity is 
very crucial. The best results are obtained within the following range of parameters: 
0.4 < L/W < 1.6, 6.25 < ^ < 16, i/W < 0.2, and h/W < 0.1. This range should 
be met naturally for most applications, and since we have based this work on creeping 
flow, the Reynolds number of the perfusion flow should be kept below or around unity. 

In all, we thus find that equation [H] can serve as a fairly accurate design guide, 
applicable for designing microfiuidic bio-reactors. 

6. Direct 3D simulation 

Up to this point we have relied on the 2D flow model based on the Brinkman- 
Darcy equation. To validate this approach and test the guideline parametrization of 
equation ([HI) , we made a full 3D direct numerical simulations of the derived bio- 
reactor design for a given set of parameters. The resulting system was modeled and 
solved in COMSOL using a ordinary laptop computer, and the solution is presented in 
figure [HI where both iso-surfaces of the pressure and streamlines are showed inside the 
computational domain. 

Similar to the earlier quasi 3D solutions of the optimized design, the pressure is 
nicely homogenized in the region of main reactor height, and also the streamlines arrange 
parallel through the main reactor. We take these results as a clear validation of the 
lubrication approach used in this work. Besides, the homogeneous flow produced by 
the design in figure [8] emphasizes the value of derived design and the parameterizations 
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Figure 8. Direct 3D numerical solution of an optimized design following the 
parametrization guide of equation 1141 Pressure iso-surfaces are gray, and streamlines 
are solid lines. Parameters are L = 0.95W,£ = O.IW, hi = 0.0414^, and A^A. 

guide of equation (fT^ . 
7. Conclusion 

To increase the utilization of continuously feed microfluidic bio-reactors, we have 
optimized the fiow-geometry of the reactor as to expose all immobilized organisms or 
substances to a very homogeneous flow field. From this we have derived a general guide- 
line of how to construct the optimal design for a broad range of reactor-dimensions. 

As the overall height of the system is much smaller than the remaining physical 
dimensions, the 3D fluid flow can essentially be described as a 2D fluid flow using a 
lubrication theory approach, where an additional volume-force arise from the viscous 
drag by the upper and lower channel- walls. 

In this work we first achieved an optimal flow-geometry by applying the free-form 
method of topology optimization. As the resulting shape in the design had a simple 
single-connected topology, we subsequently applied shape-optimization to obtain the 
different optimal geometries for various reactor-dimensions. From this analysis, we have 
constructed a general parametrization of an optimal design, which has been validated 
by direct 3D simulations. 

The design produces the homogeneous flow with a very low pressure drop, and 
this will dramatically reduce the power needed to drive the perfusion flow through the 
system. This opens the possibilities of driving the perfusion in radically new ways 
e.g. by buoyancy effects. Furthermore the fluid-volume of the flow-homogenizing design 
is minimized, which is essential when dealing with very limited fluid-samples. 

Besides applying the design to bio-reactors, it is also applicable to many other 
microfluidics system requiring perfusion of a large squared area, such as DNA 
and protein microarrays and investigation of tissue slices using fluorescent in situ 
hybridization or immuno chemistry, where samples typically are limited. 
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